function fig4_16

n=100;

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.015,'ML',0.05,'P',0.04)

ep=0.01;  w=0.5;
nn(1)=0;  nn(2)=1;
ye(1)=0;  ye(2)=1;
for i=3:n
    nn(i)=i-1;
    ye(i)=2*ye(i-1)-ye(i-2)-2*w^2*(ye(i-1)+ep*ye(i-1)^3);
end

phi=acos(1-w^2);
%th0=pi/2
%A0=1/cos(phi+th0)
aa = (ye(1) - ye(2))^2 + 2*w^2*ye(1)*ye(2);
bb = w^2*(2-w^2);
A0 = sqrt(aa/bb);
th0 = acos(ye(1)/A0);
ang = ye(1)*(1-w^2)-ye(2)
if ( ang<0 && ye(1)>0 )
    th0 = th0 -pi
elseif ((ang<0)&&(ye(1)<0))
    th0 = th0 + pi/2;
end
alpha = 0.75*w*A0^2/sqrt(2-w^2);
th0 = th0 -pi

for i=1:n
    ya(i)=A0*cos((phi+ep*alpha)*(i-1)+th0);
end

hold on
box on
grid on

plot(nn,ye,'ob','Linewidth',1.2,'MarkerSize',5)
plot(nn,ya,'r','Linewidth',1.2)
axis([0 100 -2 3.3])

set(gca,'ytick',[-2  0  2])
set(gca,'YTickLabel',{'$-2$';'0';'2'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

xlabel('$x$\,-\,axis','Interpreter','latex','FontSize',14)
ylabel('Solution','Interpreter','latex','FontSize',14)

legend(' Exact Solution',' Asymptotic Approximation','Location','North','FontSize',12)

%exportgraphics(gcf,'/Users/mark/Desktop/fig312.pdf')




